Problem Set 5 Key: Model Selection

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Tip

Try the problems yourself before consulting this key. Every number below was checked in R 4.5.1 (MuMIn 1.48.11, car 3.1.3, MASS 7.3.65) against data/wild_wheat.csv (54 genotypes, 23 columns) on 2026-10-25.

library(ggplot2)
library(MuMIn)   # model.sel(), dredge(), sw(), model.avg()
library(car)     # vif()
library(MASS)    # stepAIC()
wheat <- read.csv("https://greymonroe.github.io/PLS_206/data/wild_wheat.csv")

Part 1: Three Hypotheses

Q1. log L = −173.38, k = 3, so AIC = −2 × (−173.38) + 2 × 3 = 352.75, the same as AIC(). k counts every parameter the model estimates: the intercept, the slope of bio1, and the residual standard deviation σ (6.11 days). Students who use k = 2 get 350.75; take a point off and point them to σ.

m_temp <- lm(Flowering ~ bio1, data = wheat)
logLik(m_temp)                               # -173.3765 (df = 3)
-2 * as.numeric(logLik(m_temp)) + 2 * 3      # 352.753
AIC(m_temp)                                  # 352.753

Q2. Temperature alone is best supported: AICc 353.2, weight 0.61. Temperature + drought and temperature + day length each have weight 0.19 (ΔAICc 2.31); drought alone (ΔAICc 31.7), day length alone (38.5) and the null (50.5) have weights of essentially 0. The temperature model is about 3.2 times better supported than either runner-up (0.614 / 0.193).

The two-predictor models fit no better: their log-likelihoods (−173.364 and −173.365) are within 0.012 of the temperature model’s (−173.376). Each has one more parameter (k = 4), so AIC charges them 2 more for almost no gain in fit (AICc a little more, 2.31, because of its small-sample correction). The gap is the penalty for an extra parameter that explains nothing. Neither drought nor day length adds anything once temperature is in: their slopes are near zero with p = 0.88 (bio18) and 0.88 (Latitude). Both look important on their own only because they are correlated with temperature.

m0         <- lm(Flowering ~ 1, data = wheat)
m_dry      <- lm(Flowering ~ bio18, data = wheat)
m_day      <- lm(Flowering ~ Latitude, data = wheat)
m_temp_dry <- lm(Flowering ~ bio1 + bio18, data = wheat)
m_temp_day <- lm(Flowering ~ bio1 + Latitude, data = wheat)

tab <- model.sel(m0, m_temp, m_dry, m_day, m_temp_dry, m_temp_day)
tab
#            df   logLik  AICc delta weight
# m_temp      3 -173.376 353.2  0.00  0.614
# m_temp_dry  4 -173.364 355.5  2.31  0.193
# m_temp_day  4 -173.365 355.5  2.31  0.193
# m_dry       3 -189.219 384.9 31.68  0.000
# m_day       3 -192.644 391.8 38.54  0.000
# m0          2 -199.743 403.7 50.49  0.000
tab$weight[1] / tab$weight[2]   # 3.2

summary(m_temp_dry)   # bio18 slope -0.031, p = 0.88
summary(m_temp_day)   # Latitude slope 0.038, p = 0.88

Q3. The ranking is the same under all three. The temperature model’s weight rises from 0.57 (AIC) to 0.61 (AICc) to 0.78 (BIC), because each criterion charges more for the extra parameter in the two runner-up models: AIC charges 2 per parameter, AICc adds a small-sample correction (2k(k + 1)/(n − k − 1) = 0.48 for k = 3 and 0.82 for k = 4 with n = 54), and BIC charges log(54) = 3.99 per parameter. So under BIC the runners-up sit 3.96 behind instead of 1.97 (AIC) or 2.31 (AICc). Report AICc: n/k = 54/4 = 13.5 is well below the usual rule of thumb of 40, and AICc converges to AIC as n grows, so there is no cost to using it. (Check values for the temperature model: AICc 353.23, BIC 358.72.)

model.sel(m0, m_temp, m_dry, m_day, m_temp_dry, m_temp_day, rank = AIC)   # weights 0.573, 0.214, 0.213
model.sel(m0, m_temp, m_dry, m_day, m_temp_dry, m_temp_day, rank = BIC)   # weights 0.784, 0.108, 0.108
AIC(m_temp) + 2 * 3 * 4 / (54 - 3 - 1)   # 353.233 = AICc(m_temp)
BIC(m_temp)                              # 358.72
log(54)                                  # 3.99

Q4.

1. Both. The models are nested (temperature alone is temperature + drought with the bio18 slope set to 0). anova(): F = 0.02, p = 0.88. ΔAIC = 1.97 in favour of temperature alone. They agree: drought adds nothing.
1. AIC only. Neither model is a special case of the other, so there is no F-test. ΔAIC = 31.7 in favour of temperature.
1. Neither. AIC compares likelihoods of the same data. The log model’s likelihood is a likelihood of log(flowering days), which are numbers around 4.6 with tiny residuals, so its AIC is on a different scale and the −146 says nothing about which model is better. Both models must have the same response (and the same rows). A plus for anyone who notes the fix: put the log model back on the days scale by adding 2 Σ log(y), which gives 356.3, slightly worse than 352.8.

anova(m_temp, m_temp_dry)             # F = 0.0245, p = 0.876
AIC(m_temp_dry) - AIC(m_temp)         # 1.97
AIC(m_dry) - AIC(m_temp)              # 31.7

m_log <- lm(log(Flowering) ~ bio1, data = wheat)
AIC(m_log)                                    # -146.1
AIC(m_log) + 2 * sum(log(wheat$Flowering))    # 356.3, comparable with 352.8

Part 2: Exploring All 19 Climate Variables

Q5. Twelve rounds of pruning, dropping in order bio10 (VIF 18,585), bio5, bio11, bio19, bio16, bio17, bio2, bio1 (VIF 212), bio14, bio13, bio15 and bio9 (8.9). Six variables are left: bio3, bio4, bio6, bio8, bio12, bio18, all with VIF below 5 (largest 4.25, bio18). bio1 was not kept. VIF measures how well a variable is predicted by the other variables, not how well it predicts flowering. Annual mean temperature is close to an average of the other temperature variables (r = 0.93 with bio10, 0.85 with bio5 and bio11, 0.70 with bio6), so it was one of the most redundant columns, even though it is the best single predictor of flowering time. The rule removes redundancy blindly; it does not know which variable you care about.

full19 <- lm(reformulate(paste0("bio", 1:19), "Flowering"), data = wheat)
coef(full19)["bio7"]   # NA

vars <- setdiff(paste0("bio", 1:19), "bio7")
repeat {
  v <- vif(lm(reformulate(vars, "Flowering"), data = wheat))
  if (max(v) < 5) break
  cat("drop", names(which.max(v)), round(max(v), 1), "\n")
  vars <- setdiff(vars, names(which.max(v)))
}
vars          # "bio3" "bio4" "bio6" "bio8" "bio12" "bio18"
round(v, 2)   # 2.58 3.02 3.31 2.53 2.73 4.25
cor(wheat$bio1, wheat[, c("bio5", "bio6", "bio10", "bio11")])   # 0.85 0.70 0.93 0.85

Q6. 64 models. Top by AICc: bio3 + bio4 + bio6 (AICc 360.0, weight 0.35 in the full table, R² 0.61). Within 2 AICc units: only one other model, bio3 + bio4 + bio6 + bio12 (ΔAICc 1.43). By BIC the top model is the same, but the runner-up is the smaller bio4 + bio6 (ΔBIC 2.51), and the four-variable model drops to third (ΔBIC 2.88): BIC’s heavier penalty favours fewer predictors. Note that subset(d, delta < 2) rescales the weights to sum to 1 within the subset (0.67 and 0.33).

The predictor with the smallest p-value in the top model is bio6 (minimum temperature of the coldest month): slope −2.37 days per °C, p = 6.9 × 10⁻⁹. Genotypes from sites with milder winters flower earlier. (bio4, temperature seasonality, is close behind, p = 7.1 × 10⁻⁸; a student who plots it is fine if the band is right.)

options(na.action = "na.fail")
global <- lm(reformulate(vars, "Flowering"), data = wheat)
d <- dredge(global)
head(d)
subset(d, delta < 2)
# bio3 + bio4 + bio6          AICc 360.0  delta 0.00
# bio3 + bio4 + bio6 + bio12  AICc 361.4  delta 1.43
dredge(global, rank = "BIC")[1:3]
# bio3 + bio4 + bio6   delta 0.00
# bio4 + bio6          delta 2.51
# bio3 + bio4 + bio6 + bio12   delta 2.88

top <- get.models(d, 1)[[1]]
summary(top)   # bio3 -6.56 (p = 0.015), bio4 -0.051 (p = 7e-08), bio6 -2.37 (p = 7e-09)

p_bio6 <- ggplot(wheat, aes(bio6, Flowering)) + geom_point() + theme_classic()
p_bio6
nd <- data.frame(bio6 = seq(min(wheat$bio6), max(wheat$bio6), length.out = 100),
                 bio3 = mean(wheat$bio3), bio4 = mean(wheat$bio4))
nd <- cbind(nd, predict(top, nd, interval = "confidence"))
p_bio6 +
  geom_ribbon(data = nd, aes(bio6, ymin = lwr, ymax = upr), alpha = 0.3, inherit.aes = FALSE) +
  geom_line(data = nd, aes(bio6, fit), inherit.aes = FALSE)

Q7. Importance: bio4 1.00, bio6 1.00, bio3 0.79, bio12 0.30, bio8 0.27, bio18 0.26. A variable that explains nothing still gets a sizeable importance, because it rides along in half of the models, including copies of the good ones: the models with one useless extra term sit about 2 AICc behind and keep some weight. Here the three weakest variables cluster at 0.26 to 0.30, which is that floor, so bio12’s 0.30 is no evidence that it matters. Only bio4 and bio6 stand clearly above it, with bio3 in between.

In the average of the two models within 2 AICc, bio12’s slope is −0.073 in the conditional average (averaged only over the one model that contains it) but −0.024 in the full average. The full average counts a slope of 0 for the model without bio12 (weight 0.67), so it shrinks the slope towards 0 in proportion to how little support the variable has: −0.073 × 0.33 = −0.024. Neither average is distinguishable from zero (p = 0.66 full, 0.32 conditional). bio3, bio4 and bio6 are in both models, so their two averages are identical.

sw(d)
#   bio4 bio6 bio3 bio12 bio8 bio18
#   1.00 1.00 0.79 0.30  0.27 0.26   (each in 32 of 64 models)

ma <- model.avg(d, subset = delta < 2)
summary(ma)
# full average:        bio12 -0.0240 (p = 0.66)
# conditional average: bio12 -0.0728 (p = 0.32)

Q8. The Part 1 temperature model is better supported: AICc 353.2 vs 360.0, ΔAICc = 6.7, about 29 times the support (exp(6.7/2)). The top dredge model is only the best of the models dredge() was allowed to try, and the VIF pruning had already removed bio1, the strongest single predictor. So the dredge winner is a hypothesis (winter cold and temperature seasonality as cues for flowering), not a finding: it came from searching 64 models on these same 54 genotypes and would need new data to test. The dredge serves exploration; the Part 1 comparison of a few hypotheses stated in advance is inference.

AICc(m_temp)                      # 353.2
AICc(top)                         # 360.0
exp((AICc(top) - AICc(m_temp)) / 2)   # 29

Part 3: The Traps

Q9. (a) No. Backward from the 18-variable model stops at seven variables, bio1 + bio3 + bio11 + bio13 + bio16 + bio17 + bio18 (AIC 354.2); forward from the null stops at bio1 alone (AIC 352.8). Backward ends at a worse AIC than forward: each step only looks one variable ahead, so it can get stuck, and the answer depends on where you start. The seven-variable model is also a warning sign: most of its p-values sit between 0.02 and 0.05, the pattern you get when a search keeps whatever happens to help a little.

Kept in 200 resamples (set.seed(206)): bio6 99.5%, bio4 99%, bio3 79.5%, bio12 28.5%, bio8 19%, bio18 10.5%. bio4 and bio6 are stable; the rest come and go. The 200 resamples gave 15 different final models; the most common, bio3 + bio4 + bio6 (the Q6 winner), came up 94 times (47%), so in more than half of the resamples stepwise picked something else. The number of variables kept ranged from 2 to 6. Exact percentages depend on the seed and the loop; anything close to this pattern is correct.

full18 <- lm(reformulate(setdiff(paste0("bio", 1:19), "bio7"), "Flowering"), data = wheat)
s_back <- stepAIC(full18, trace = 0)
formula(s_back)   # bio1 + bio3 + bio11 + bio13 + bio16 + bio17 + bio18
summary(s_back)
AIC(s_back)       # 354.2
s_fwd <- stepAIC(lm(Flowering ~ 1, data = wheat),
                 scope = list(lower = ~ 1, upper = formula(full18)),
                 direction = "forward", trace = 0)
formula(s_fwd)    # Flowering ~ bio1
AIC(s_fwd)        # 352.8

set.seed(206)
kept <- matrix(0, 200, length(vars), dimnames = list(NULL, vars))
mods <- character(200)
for (b in 1:200) {
  wb <- wheat[sample(nrow(wheat), replace = TRUE), ]
  s <- stepAIC(lm(reformulate(vars, "Flowering"), data = wb), trace = 0)
  k <- attr(terms(s), "term.labels")
  kept[b, k] <- 1
  mods[b] <- paste(k, collapse = " + ")
}
colMeans(kept)          # bio3 0.795, bio4 0.990, bio6 0.995, bio8 0.190, bio12 0.285, bio18 0.105
length(unique(mods))    # 15
sort(table(mods), decreasing = TRUE)[1:3]   # bio3 + bio4 + bio6: 94
table(rowSums(kept))    # 2:12, 3:118, 4:57, 5:12, 6:1

Q10. (a) Yes. On all 54 rows, stepAIC() keeps bio3, bio4, bio6 and four noise columns, and noise9 comes out with slope −2.46, p = 0.0046. From (a) alone you would have reported that noise9 predicts flowering time, a “significant” effect of a column of random numbers.

Selecting on half A, the terms with p < 0.01 are bio4, bio6 and noise5 (p = 0.0048 in half A; noise9 did not make the cut this time, another sign that it was chance). Refit on the other 27 genotypes: bio4 p = 0.0021 and bio6 p < 0.0001, same signs, so both are confirmed; noise5 p = 0.82, gone. A p-value from a model chosen on the same data is not a real test: the search has already kept whatever happened to fit, so the p-values are too small. The honest test is on data the search never saw. This is not bad luck with one seed: across 100 noise datasets, backward stepAIC() kept at least one noise column with p < 0.05 in 52% of them (13% with p < 0.01).

A good answer may also note that some genotypes share a collection site and so have identical climate values (only 31 distinct climates among the 54 genotypes), so the two halves are not fully independent; a stricter split would put whole sites in one half.

set.seed(19)
wn <- wheat[, c("Flowering", "bio3", "bio4", "bio6", "bio8", "bio12", "bio18")]
for (j in 1:10) wn[[paste0("noise", j)]] <- rnorm(nrow(wn))
half <- sample(nrow(wn), 27)

s_all <- stepAIC(lm(Flowering ~ ., data = wn), trace = 0)
summary(s_all)   # bio3, bio4, bio6, noise1, noise5, noise7, noise9; noise9 p = 0.0046

s_A <- stepAIC(lm(Flowering ~ ., data = wn[half, ]), trace = 0)
cA <- summary(s_A)$coefficients
sigA <- rownames(cA)[-1][cA[-1, 4] < 0.01]
sigA             # "bio4" "bio6" "noise5"
m_B <- lm(reformulate(sigA, "Flowering"), data = wn[-half, ])
summary(m_B)     # bio4 p = 0.0021, bio6 p < 0.0001, noise5 p = 0.82

# Context (not asked): how often does noise end up "significant"?
set.seed(206)
hit <- replicate(100, {
  d2 <- wheat[, c("Flowering", vars)]
  for (j in 1:10) d2[[paste0("noise", j)]] <- rnorm(nrow(d2))
  cf <- summary(stepAIC(lm(Flowering ~ ., data = d2), trace = 0))$coefficients
  any(cf[grep("noise", rownames(cf)), 4] < 0.05)
})
mean(hit)   # 0.52
nrow(unique(wheat[, paste0("bio", 1:19)]))   # 31 distinct climates